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Abstract 

Isogeometric  analysis  is  applied  to  boundary  integral  equations  corresponding  to  boundary-value 
problems  governed  by  Laplace’s  equation.  It  is  shown  that  the  smoothness  of  geometric  parametriza- 
tions  central  to  computer-aided  design  can  be  exploited  for  regularizing  integral  operators.  As  a  result, 
one  obtains  high-order  collocation  methods  based  on  superior  approximation  and  numerical  integration 
schemes  and  well-conditioned  systems  of  linear  algebraic  equations.  It  is  demonstrated  that  the  proposed 
approach  allows  one  to  solve  boundary-value  problems  with  an  accuracy  close  to  machine  precision. 


1  Introduction 

Isogeometric  analysis  (IgA)  [24,  35]  is  a  framework  for  numerical  schemes  for  solving  boundary-value  problems 
(BVPs)  in  which  the  basis  functions  coincide  with  those  used  for  geometric  parametrizations  in  computer 
aided  design  (CAD).  Thus,  in  contrast  to  conventional  finite  element  methods,  IgA  relies  on  Non-Uniform 
Rational  B-splines  (NURBS)  [43,  46],  T-splines  [55,  53]  or  subdivision  surfaces  [21,  48,  51]  rather  than  piece- 
wise  polynomials  as  the  basis  functions.  IgA  simplifies  mesh  generation,  and  thus  it  significantly  shortens 
the  design-through-analysis  process  for  high-end  engineering  components.  Furthermore,  IgA  is  advantageous 
because,  in  contrast  to  finite  element  methods,  it  fully  preserves  geometry  of  CAD-generated  shapes  and 
involves  basis  functions  with  attractive  properties.  These  features  have  given  rise  to  accurate  and  efficient 
numerical  schemes  successfully  applied  to  fluid  mechanics  [1,  3,  5,  7,  8,  12,  29],  solid  mechanics  [39],  elec¬ 
tromagnetism  [20],  fluid-structure  interaction  [9,  10,  11,  60],  structural  dynamics  [25,  26],  plates  and  shells 
[15,  16,  27,  28,  37,  22,  23],  phase-field  models  [17,  32,  33],  and  shape  optimization  [40,  41,  45,  59]. 


Boundary  element  methods  (BEMs)  are  numerical  schemes  for  solving  boundary  integral  equations  (BIEs). 
Like  finite  element  methods,  BEMs  rely  on  piecewise  polynomials  for  approximating  the  geometry  and  field 
variables.  Thus,  by  replacing  piecewise  polynomials  with  NURBS  or  T-splines,  one  can  develop  isogeometric 
BEMs.  This  approach  has  been  already  undertaken  [13,  14,  31,  38,  42,  44,  54,  56]. 


The  premise  of  this  paper  is  that  IgA  can  radically  improve  numerical  schemes  for  solving  BIEs  because  of 
the  additional  smoothness  of  NURBS  and  T-splines  in  comparison  to  C°-continuous  piecewise  polynomials. 
We  show  that  one  can  regularize  the  key  singular  integral  operators,  and  construct  superior  approximation 
and  integration  schemes  and  well-conditioned  systems  of  linear  algebraic  equations.  These  schemes  allow 
one  to  solve  BIEs  with  machine  precision.  We  demonstrate  these  advantages  of  IgA  by  applying  it  to 
BIEs  corresponding  to  Laplace’s  equation.  This  restriction  significantly  simplifies  mathematical  aspects  of 
our  work,  but  most  of  the  results  can  be  readily  extended  to  equations  of  classical  elasticity  and  Stokes 
equations  of  fluid  mechanics. 
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In  the  context  of  BEMs,  the  proposed  approach  can  be  classified  as  a  high-order  collocation  method  allowing 
for  weakly-singular,  singular,  and  hyper-singular  integral  operators.  Previously  this  setting  was  possible 
with  Galerkin  but  not  collocation  methods.  Accordingly,  we  exploited  analytical  tools  different  from  those 
usually  used  in  analysis  of  Galerkin  methods  [49,  58].  In  this  regard,  we  believe  that  our  approach  may  be 
useful  for  developing  mathematical  foundations  for  collocation  schemes  for  BIEs.  In  addition,  collocation 
BEMs  dominate  engineering  applications,  thus  we  believe  that  our  results  are  not  only  of  mathematical  but 
also  practical  interest. 


The  rest  of  the  paper  is  organized  as  follows.  In  Section  2,  we  define  the  model  BVP,  corresponding  BIEs, 
and  a  proper  continuous  setting  for  BVPs  defined  on  domains  with  smooth  boundaries.  In  Section  3, 
we  introduce  collocation  schemes  for  the  integral  operators,  defined  on  smooth  surfaces.  In  Section  4, 
we  demonstrate  that  all  collocated  integral  operators  defined  on  smooth  surfaces  can  be  reduced  to  weakly- 
singular  integrals.  In  Section  5,  we  briefly  summarize  basic  results  from  IgA  and  identify  extensions  needed  for 
solving  BIEs.  We  also  describe  the  extension  of  our  methodology  to  piecewise  smooth  surfaces  which  allows 
us  to  accommodate  multiple  patches  with  (^-continuity,  degenerate  patches  resulting  in  local  C-continuity, 
and  extraordinary  points  or  star-points.  In  Section  6,  we  consider  representative  example  problems,  which 
allow  us  to  demonstrate  various  important  mathematical  and  computational  aspects.  In  Section  7,  we 
summarize  key  results  of  this  work  and  briefly  discuss  directions  for  future  research.  Mathematical  details 
are  presented  in  the  appendix. 


2  Continuous  Formulation 

2.1  Model  boundary- value  problem 

Consider  a  bounded  domain  O  C  JR3  with  T  :=  <90.  We  assume  that  T  is  a  C2-surface  (r  £  C2),  which, 
loosely,  means  that  P  can  be  mapped  on  R2,  and  the  inverse  of  that  map  ip  is  twice  continuously  differentiable. 
For  a  rigorous  definition  of  C2-surfaces  we  refer  to  the  appendix.  This  restriction  T  £  C2  is  sufficient  for 
establishing  mathematical  foundations  for  integral  equations.  A  typical  CAD  parametrization  map  can  be 
defined  as  a  union  of  C2 -surfaces,  resulting  in  a  globally  Lipschitz  T;  we  denote  this  class  of  surfaces  by 
C2 .  Unfortunately,  mathematical  foundations  for  BIEs  defined  on  C2-surfaces  are  not  well  developed  at  this 
stage.  Thus,  in  this  and  two  following  sections,  we  restrict  our  attention  to  T  £  C2.  On  the  other  hand,  the 
loss  of  smoothness  of  C-surfaces  can  be  compensated  by  using  appropriate  numerical  schemes  presented  in 
Sections  5  and  6. 


The  model  BVP  is  formulated  for  Laplace’s  equation 

—A u  =  0  in  D,  (1) 

and  mixed  boundary  conditions,  which  include  Dirichlet 

u  =  gD  on  Td,  (2) 


and  Neumann 


t:=n-'S7u  =  gN  on  Tjv  (3) 

data.  Here  n  denotes  the  outward  unit  normal  vector  onT;  T^UTjv  =  T  and  Td  (~l  Tn  =  0.  At  this  stage, 
we  require  T  £  C2,  but  we  do  not  specify  the  smoothness  of  u,  go,  and  g n.  This  will  be  done  once  we 
introduce  the  integral  operators. 
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2.2  Integral  equations 

The  integral  equations  equivalent  to  (l)-(3)  and  their  mathematical  properties  are  well  known  [34,  49,  58]. 
They  involve  the  representation  formula  that  allows  one  to  determine  the  solution  u  in  terms  of  the  Cauchy 
boundary  data  (u,t): 


u{x)  =  j  G(x,  y)t(y)dsy  -  j  n(y )  •  [VyG(x,  y )]  u(y)dsy  x  €  H  , 


(4) 


where 


G(x,y)  = 


47r|a!  —  y\ 


is  the  fundamental  solution  of  Laplace’s  equation.  The  Cauchy  data  can  be  reconstructed  using  the  singular 
boundary  integral  equation  (SBIE) 


Here  X  is  the  identity, 


is  the  single-layer  operator,  and 


^X  +  K,  )  u  =  Vt  on  T  . 


Vt(x)  :=  J  G(x,y)t.(y)dsy,  xeT 


lCu(x)  :=  J  n(y)  ■  [VyG(x,  y)}  u(y)dsy,  i£f 


is  the  double-layer  operator.  Alternatively,  the  Cauchy  data  can  be  reconstructed  using  the  hyper-singular 
integral  equation  (HSBIE) 

1 


-I  —  K,'  )  t  =  Vu  on  T 


where 


K't{x)\=  J  n(x)  ■  [VxG(x,y)\t(y)dsy,  xeT 
is  the  adjoint  double-layer  operator,  and 

Vu(x)  :=  -  J  n(x)  ■  {Vxn(y)  ■  [\7yG{x,  y)}}  u(y)dsy,  x  G  T 
is  the  hyper-singular  operator. 


2.3  Analysis  of  integral  equations 

Since  T  £  C2  one  can  prove  that  the  operators 

1C, 1C1  :  C(r)  C(T) 

are  compact;  see  Appendix.  Here  C'(T)  is  the  space  of  continuous  functions  on  T.  Similarly,  we  adopt 
the  same  convention  for  other  spaces  such  as  C2(r)  which  is  the  space  of  twice  continuously  differentiable 
functions  on  T.  As  a  consequence  of  the  compactness,  Fredholm’s  alternative  implies  that  the  operators 

Qt  +  k)  :  C*(r)  -»•  C«(r),  and  Qj  -  K^j  :  C(T)  ^  C'(T) 
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are  invertible.  Here  C*(T)  is  the  space  of  all  functions  u  G  C(T)  with 

J  u(y)dsv  =  0  , 

and  C**(P)  is  the  space  of  all  functions  u  G  C(T)  with 

J  u{y)weq{y)dsv  =  0  , 

where  weq  =  V_1l  is  the  natural  weight,  which  has  been  used  in  analysis  of  the  hyper-singular  operator  [58]. 
Further,  compactness  of  K,  and  1C'  implies  that  SBIE  and  HSBIE  should  be  solved  by  inverting  \X  +  K,  and 
\X  —  K.' ,  respectively. 


For  the  pure  Neumann  BVP,  the  SBIE  takes  the  form 

'1 


X  +  K,  )  u  =  VgN  on  T. 


(5) 

Since  V  :  C(P)  — ►  C(r)  and  the  boundary  data  g n  must  satisfy  the  solvability  condition 

J  gN(x)dsx  =  0  , 

it  follows  that  g n  G  C(r)  implies  VgN  G  C'**(r).  Therefore,  the  invertibility  of  +  K.)  implies  that  (5) 
has  a  unique  solution  for  u  G  C*(F). 


For  the  pure  Diriclilet  BVP,  the  HSBIE  takes  the  form 

(^X-IC^t  =  VgD  onF.  (6) 

One  can  prove  that  go  G  C2(r)  implies  Vgo  G  0(F).  Thus  the  invertibility  of  (\X  —  K.')  implies  that  (6) 
has  a  unique  solution  for  t  G  C(r). 


Based  on  analysis  of  the  pure  BVPs,  it  is  appropriate  to  require  go  G  C2{Td)  and  g n  G  C(Pjv)  for  the 
mixed  BVP.  These  restrictions,  however,  do  not  guarantee  u  G  C2(r)  and  t  G  O(r).  To  formulate  the  BIEs 
corresponding  to  the  mixed  BVP,  we  define  extensions  go  G  C2(r)  and  g n  G  O(r),  and  express  the  Cauchy 
data  as 


By  construction,  «.|rD 
and  t: 


and 


0  and  t|rN 


u  =  u  +  go  and  t  =  t  +  gN  ■  (7) 

=  0.  Now  we  can  rewrite  the  SBIE  and  HSBIE  as  the  equations  for  u 


Vi  -  (  \xi  +  ic  \  u  =  (^I  +  IC)gD-  VgN  on  T, 


(8) 


t-Vu  =  VgD  -  (^X  -  K,'  )  gN  on  T 


(9) 


Note  that  while  gp  and  gpf  are  not  uniquely  defined,  the  structure  of  (8)  and  (9)  is  such  that  u  and  t  are 
uniquely  defined  as  long  as  (8)  and  (9)  are  uniquely  solvable  for  u  and  t.  Since  go  G  C2(P)  and  g n  G  C'(r), 
the  right-hand  sides  of  both  equations  are  in  C(P).  However,  this  is  insufficient  for  establishing  unique 
solvability  for  (8)  and  (9),  under  the  provisions  u  G  C2(r)  and  t  G  C'(P). 
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Remark  2.1.  It  is  straightforward  to  extend  our  analysis  to  the  pure  Robin  BVP  in  which  the  boundary 
data  are  prescribed  as 


t  +  ku  =  ga  on  r  , 

where  k  €  L°°(r)  and  g r  is  a  prescribed  function  in  C(T).  This  problem  is  similar  to  the  pure  Neumann 
problem. 


Remark  2.2.  The  regularity  requirements  on  T,  go  and  gw  can  be  slightly  relaxed.  The  results  of  this  section 
can  be  extended  to  continuously  differentiable  r,  gjy  with  Lipschitz  continuous  derivatives,  and  gw  €  L°° . 

Remark  2.3.  While  it  appears  that  the  condition  go  €  C2(T)  is  too  restrictive,  in  practice  go  is  often  a 
constant.  For  example,  in  the  context  of  heat  conduction,  Dirichlet  boundary  conditions  represent  a  situation 
in  which  is  placed  in  a  constant-temperature  environment  whose  temperature  is  not  affected  by  H.  If  it  is 
the  case,  a  constant  go  is  simply  extended  to  the  entire  boundary,  so  that  go  is  constant. 

3  Collocation  discretization 

Let  us  consider  the  mixed  BYP.  Approximations  for  u  and  t  are  constructed  as 

D  N 

n  n 

Uh(x)  =  ^2  &[A]NA  (x)>  and  ih{x)  =  t[A]NA  ix)  1 
.1-1  j4= 1 

respectively,  where  NA  (x)  and  Nf  (x)  are  the  basis  functions  and  u  and  t  are  column- vectors.  Since  u(x)  =  0 
for  x  G  Td  and  t(x)  =  0  for  x  €  Tjv,  the  basis  functions  are  such  that  NA(x)  =  0  for  x  €  To  and  NA(x)  =  0 
for  x  £  T]sf.  Accordingly,  we  define  collocation  points  xA  on  and  x A  on  This  assignment  of  the 
superscripts  may  be  somewhat  confusing,  but  it  simply  reflects  the  fact  that  the  basis  functions  NA[x)  and 
NA  (x)  are  supported  on  Tn  and  T ,  respectively. 


Upon  collocating  (8)  at  xA  and  (9)  at  xA,  one  generates  the  system  of  linear  algebraic  equations  for  u  and 
t: 

Vt-(J2ID  +  K^u  =  ls  (10) 

and 

(^IN-K^f-Du  =  lH.  (11) 

Here  the  components  of  the  system  matrices  are  defined  as 


V[A,  B]  :=  VN» (: r E)  =  /  G(x%,y)N% {y)dsv  ,  (12) 

Jtd 

Id[A,B]:=NE(x%),  (13) 

I<[A,B]  :=  ICNE(x^)  =  [  nfy)  •  [VyG(xE,y)\NE  (y)dsy  ,  (14) 

In[A,B]:=NE(x%),  (15) 

K'[A,B]  :=1CNE(x%)=  f  n(xNA)  •  [V  xG(x%,y)}N»  (y)dsy  ,  (16) 

JrD 

D[A,B]  :=VNE(x%)  =  -  [  n  (a£)  •  {V*  [n(y)  ■  VyG(x%,  y)\  }  NE(y)day  .  (17) 

JrN 
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The  right-hand-side  vectors  are  defined  as 


fs[A]  ■=  ^9d{xa)  +  ^9d{xa)  -  VgN( x%)  , 

Lh[A]  ■=  v9d(x^)  -  ^gN(x%)+lCgN{x%) . 


For  the  pure  Neumann  BVP,  all  collocation  points  are  x%,  and  therefore  it  is  sufficient  to  solve  (10)  by 
setting  th  =  0,  gD  =  0,  and  gN  =  gN- 

-{lID+K)^  =  ls-  (18) 

The  compactness  of  K,  can  be  used  to  show  that  (18)  is  uniquely  solvable  as  long  as 

J  Uh(x)dsx  =  0  .  (19) 

Furthermore,  under  condition  (19),  Uh  converges  to  the  exact  solution  u  of  (5)  at  optimal  rates. 


For  the  pure  Dirichlet  BVP,  all  collocation  points  are  ar^ ,  and  therefore  it  is  sufficient  to  solve  (11)  by  setting 
uh  =  0,  gN  =  0,  and  gD  =  gD: 


{]2lN-K')t  =  lH.  (20) 

The  compactness  of  1C  can  be  used  to  show  that  (20)  is  uniquely  solvable  and  th  converges  to  the  exact 
solution  t  of  (6)  at  optimal  rates. 


Remark  3.1.  It  is  straightforward  to  extend  the  collocation  scheme  to  the  pure  Robin  BVP  following  the 
prescription  for  the  Neumann  BVP.  Also,  for  the  pure  Robin  BVP,  unique  solvability  and  optimal  convergence 
rates  can  be  established  in  a  way  similar  to  the  pure  Neumann  B  VP. 

Remark  3.2.  In  all  cases,  the  system  of  governing  algebraic  equations  is  constructed  so  that  one  has  to  invert 
matrices  associated  with  the  operators  |X+/C  and  \T—KI.  This  construction  results  in  well-conditioned  linear 
algebraic  systems  [2],  and  it  is  superior  to  alternative  formulations;  see  Section  6  for  numerical  examples. 

4  Regularization  of  operators 

In  general,  the  SBIE  and  HSBIE  involve  integrals  containing  weakly-singular,  singular,  and  hyper-singular 
kernels.  In  this  section,  we  establish  that  for  T  £  C2  and  sufficiently  smooth  functions,  as  established  in 
Section  2,  all  integrals  can  be  evaluated  as  weakly  singular  ones.  This  property  is  essential  as  it  allows  one 
to  develop  numerical  integration  schemes  with  spectral  accuracy  [52]. 


The  single-layer  operator  V  is  naturally  weakly  singular  for  the  continuous  data;  the  same  is  true  for  K,  and 
1C  (see  Appendix).  To  regularize  the  hyper-singular  operator  T>  we  begin  with  approximating  u  £  C2(P)  in 
the  vicinity  of  x: 

u{y)  =  u(x)  +  (VTu)(x)  -(y-x)  +  0( \x  -  y\ 2) ,  (21) 
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where  (\7tu)(x)  is  the  tangential  gradient  of  u  on  T.  With  this  approximation,  the  hyper-singular  operator 
can  be  expressed  as 


Vu(x)  =  -  j  n(x)  ■  X7x  {n(y)  ■  [VyG(x,  y)]}  u(y)dsv 

=  -  j  n{x)  ■  \7X  {n(y)  ■  [VyG{x,  y)]}  [ u(y )  -  u{x)  -  ( VTu){x )  •  (y  -  x)\dsy 
-  u(x)  J  n{x)  ■  Vx  { n{y )  ■  [V yG(x,  y)]}  dsy  -  J  n{x)  ■  Wx  { n(y )  •  [VyG(x,  y)]}  (VTu)(x)  ■  (y  -  x)ds 


y  • 
(22) 


Since  constant  and  linear  functions  are  harmonic,  the  HSBIE  implies  that 
-  J  n(x )  •  Vx  (n(y)  •  [VyG(x,  y)]}  dsy  =  0, 

-  J  n(x)  ■  {V xn(y)  ■  [VyG( x,  y)]}  (VTu)(x)  ■  (y  -  x)dsv  =  -  j  n(x)  ■  V x  [G(x,  y)]  n(y)  ■  (VTu)(x)dsy. 


Now  the  last  line  of  (22)  can  be  rewritten  as 


Vu(x) 


n(x)  ■  Vx  {n(y)  •  [VyG(x,  y)]}  [u(y) 


u(x)  -  (\/Tu)(x)  ■  (y  -  x)\dsy 


n(x)  ■  S7X  [G(x,y)\  n(y)  ■  ( VTu)(x)dsy  . 


In  this  equation,  the  first  term  on  the  right-hand  side  is  weakly  singular  because  of  (21)  and  the  second  term 
is  weakly  singular  because  it  is  equal  to  JC[n(y)  ■  (\7Tu)(x)}.  Our  development  closely  follows  that  in  [36]. 
An  important  aspect  of  the  regularization  scheme  is  computing  the  tangential  gradient.  This  issue  will  be 
addressed  in  Section  5  once  IgA  parametrizations  have  been  introduced. 


5  Isogeometric  Analysis 

5.1  CAD  geometry 

In  IgA  it  is  presumed  that  the  surface  T  is  described  using  a  CAD  tool.  Invariably  those  descriptions  rely 
on  B-splines.  A  one-dimensional  B-spline  of  degree  p  is  a  piecewise  polynomial  function  of  degree  p.  The 
smoothness  between  polynomials  can  be  controlled  locally  and  can  be  up  to  Gp_1.  Thus,  if  desired,  one 
can  construct  a  B-spline  of  degree  p ,  which  is  Gp_1  globally.  Multi-dimensional  B-splines  are  constructed  as 
tensor  products  of  one-dimensional  B-splines.  T-splines  are  constructed  as  linear  combinations  of  B-splines 
without  the  need  for  the  tensor-product  structure.  This  feature  is  very  attractive,  as  it  allows  local  refinement 
with  hanging  nodes,  and  therefore  we  use  T-splines.  For  further  details  on  T-splines  we  refer  to  [53,  54], 


First,  let  us  assume  that  T  can  be  mapped  on  a  rectangular  parametric  domain  F  C  R2.  In  CAD  the  map 
(f  :  r  — >  r  is  constructed  in  terms  of  a  set  of  control  points  Pa  €  R3,  weights  wa  >  0,  and  T-splines  AJ 
defined  on  F: 


x  =  <A£i,6) 


EaPawaN^^) 

1A2) 


(6,6)  ef. 


(23) 


This  map  has  numerous  advantages  over  the  simpler  map  6)-  In  particular,  it  can  represent 

quadric  surfaces  exactly.  Further,  since  T-splines  AJ  are  a  superset  of  B-splines,  the  map  (23)  can  be 
restricted  to  NURBS,  which  are  currently  an  industrial  standard.  Note  that  in  principle  the  map  smoothness 
can  be  controlled  by  choosing  appropriately  smooth  T-splines.  However,  in  practice  it  is  standard  to  set 
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p  =  3  and  use  C2  T-splines,  so  that  ip  £  C2(f).  By  adopting  a  global  definition  for  the  map  ^  :=  <p,  we 
conclude  that  T  defined  by  p  is  a  C2-surface. 


The  map  defined  in  (23)  requires  a  rectangular  f  and  therefore  it  is  rather  limited.  For  example,  it  cannot 
be  used  for  constructing  a  cube.  In  the  context  of  NURBS,  this  issue  is  usually  addressed  by  allowing  the 
parametric  domain  to  consist  of  multiple  rectangular  patches.  This  creates  a  new  host  of  problems  associated 
with  imposing  continuity  conditions  across  patches.  This  issue  is  naturally  resolved  with  T-splines,  as  they 
allow  for  hanging  nodes  and  smooth  basis  functions  across  patches.  Nevertheless,  even  with  T-splines,  one 
has  to  address  extraordinary  points.  By  definition,  those  points  are  intersections  of  three  or  more  than  four 
patches.  At  extraordinary  points,  the  parametrization  map  p  is  only  C°,  and  as  a  result  T  £  C2.  For 
details  we  refer  to  [54] .  Another  way  of  generalizing  (23)  is  by  allowing  rectangles  to  be  mapped  on  triangles 
by  collapsing  edges.  This  approach  involves  two  ingredients:  (i)  in  the  parametric  domain,  all  T-splines, 
supported  on  the  edge  to  be  collapsed,  are  constructed  as  C°  functions  across  the  edge,  and  (ii)  control  points 
corresponding  to  T-splines  supported  on  the  edge  to  be  collapsed  are  assigned  to  the  same  position.  Like  the 
treatment  of  extraordinary  points,  this  construction  yields  locally  C'°-parametrizations.  Thus  generalized 
maps,  involving  multiple  rectangular  patches,  give  rise  to  F  £  C2. 


parametric  space  physical  space 


f  r 


Figure  1:  Parametric  and  physical  spaces  for  a  torus. 


5.2  Basis  functions 

Let  us  suppose  that  control  points  Pa,  weights  wa,  and  T-splines  iVj,  prescribing  T  via  (23)  are  given.  Then 
in  the  parametric  domain  the  basis  functions  are  constructed  using  the  partition  of  unity 


NA(t  i,6) 


i,g2) 


(24) 


The  basis  functions  Na  in  the  physical  domain  T  are  constructed  via  the  standard  map 

Na  :=  Na  o  p^1  . 


(25) 


This  construction  includes  extraordinary  points  but  not  collapsed  edges.  For  the  latter  cases,  the  basis 
functions  Na'  supported  on  a  collapsed  edge  are  replaced  by  a  single  basis  function 

n  =  J2Na>- 

A' 
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After  that,  the  corresponding  basis  function  N  is  constructed  via  (25).  Note  that  in  our  work,  both  extraor¬ 
dinary  points  and  collapsed  edges  give  rise  to  basis  functions  which  are  locally  C 2  but  globally  continuous. 
For  extraordinary  points,  one  could  use  a  constrained  optimization  framework  that  gives  rise  to  C^-basis 
functions;  for  details  see  [54]. 


For  solving  BIEs,  one  also  needs  discontinuous  basis  functions  for  approximating  t,  which  can  be  discon¬ 
tinuous  and  even  singular  even  if  T  £  C2  and  prescribed  Cauchy  data  are  smooth.  It  is  straightforward  to 
define  discontinuous  T-splines  iVj’dlsc.  However,  CAD  parametrizations  involve  the  weights  for  continuous 
T-splines  only.  While,  in  principle,  one  can  compute  the  weights  for  discontinuous  T-splines,  and  then  use 
(24),  we  adopt  a  simpler  construction  involving  unweighted  scaled  basis  functions: 


^a(£i,£2) 


6) 

£b=i  wbN%(Zi,&)  ' 


(26) 


The  scaling  by  y 
struction  in  (26) 


n  -i 


£b=i^I(5i,6)J  is  motivated  by  numerical  examples  rather  than  theory.  The  con- 
is  not  a  partition  of  unity,  but  this  property  is  not  required  for  analysis  of  BIEs. 


5.3  Collocation  points 

The  Greville  abscissa  of  a  B-spline  is  a  point  in  the  parametric  domain  T  whose  coordinates  are  defined  as 
the  average  of  the  coordinates  of  the  knots.  These  points  often  correspond  to  the  maximum  value  of  the 
B-spline.  It  has  been  shown  that  they  are  ideally  suited  for  interpolation  and  that  they  can  be  naturally 
extended  to  T-splines.  Further,  Greville  abscissae  have  been  widely  used  as  collocation  points  for  the  basis 
functions  (25)  in  various  numerical  methods  [6,  4,  50],  including  those  for  BIEs  [38,  54,  56,  57].  However,  it 
has  been  recognized  [54]  that  for  discontinuous  T-splines  Greville’s  abscissae  may  coincide,  and  therefore  one 
needs  to  modify  the  construction.  This  issue  has  been  addressed  by  introducing  “2-ring”  collocation  points 
[54]  for  discontinuous  cubic  T-splines. 


In  this  work,  we  generalize  the  definition  of  the  “2-ring”  collocation  points  to  T-splines  of  degree  p.  To  this 
end,  let  us  consider  a  one-dimensional  B-spline  H(£)  of  degree  p  with  a  <  £  <  b.  If  B(£)  £  C  ([a,  6]),  then 
the  “2-ring”  collocation  point  is  simply  Greville  abscissa.  If  £?(£)  is  discontinuous  at  £  =  a  then  the  “2-ring” 
collocation  point  is 


£2 


—ring 


:=  a 


b  —  a 

p  +  2  ' 


if  the  discontinuity  is  at  £  =  b ,  then 


£2 


—ring 


:=b  + 


a  —  b 

p  +  2  ' 


To  construct  the  collocation  point  for  a  two-dimensional  discontinuous  T-spline,  we  can  exploit  that  locally 
(rather  than  globally)  T-splines  are  tensor  products  of  one-dimensional  B-splines.  Therefore  once  a  two- 
dimensional  discontinuous  T-spline  is  represented  by  Tij(£  1,^2)  =  £?i(£i  )-£!,•  (£2),  one  can  find  the  coordinates 
of  the  “2-ring”  collocation  point  by  treating  B,:(£i)  and  Hj(£ 2)  separately. 


5.4  Numerical  Integration 

It  has  been  established  in  Section  4  that,  upon  collocation,  all  operators  can  be  evaluated  as  weakly-singular 
integrals  on  C2-surfaces.  In  this  section,  we  focus  on  numerical  integration  schemes  applicable  to  the  opera¬ 
tors  on  Cl2-surfaces.  First,  let  us  establish  that  all  operators  can  be  evaluated  as  weakly-singular  integrals  on 
C2-surfaces,  as  long  as  approximations  are  allowed  to  include  discontinuous  basis  functions.  This  provision 
is  essential  as  it  allows  one  to  move  the  collocation  points  away  from  surface  irregularities,  associated  with 
either  CAD  parametrizations  of  smooth  surfaces  (extraordinary  points  and  collapsed  edges)  or  non-smooth 
surfaces. 
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Figure  2:  Greville  abscissae  (large  grey  circles)  and  “2-ring”  Greville  abscissae  (small  black  circles)  for  a 
single  patch  in  (a)  1-D  and  (b)  2-D.  Note  that  the  abscissae  coincide  except  for  the  patch  boundaries. 


For  discontinuous  basis  functions  the  collocation  points  are  restricted  to  the  smooth  part  of  T.  Then  among 
the  operators  defined  in  (12),  K!  and  V  can  be  evaluated  using  weakly-singular  integrals  because  they  operate 
on  functions  that  are  smooth  in  neighborhoods  of  the  collocation  points.  Then  local  Taylor  expansions  can 
be  exploited  for  regularization,  similar  to  the  way  it  is  done  in  Section  4.  Further,  the  operator  V  is  naturally 
weakly  singular  on  C'2-surfaces.  To  establish  weak  singularity  of  (al  +  K.)  u ,  let  us  substitute  u(x)  =  1  in 
the  SBIE  to  obtain 

(crX  +  JC)  1  =  0 

and  the  regularization 

(crZ  +  K)  u(x)  =  j  n(y)VyG(x,  y)u{y)  [u{y)  -  u(x)\  dsy  . 

Note  that  we  replaced  1/2  with  a  to  reflect  the  fact  that  F  is  not  smooth  near  certain  x.  This  regularization 
is  not  sufficient  for  making  the  integral  weakly  singular  for  u  £  C(r).  On  the  other  hand  if  u  is  Lipschitz 
continuous,  so  that 

I u(y)  -  u(x)  |  <C\y-x\  , 

then  the  integral  becomes  weakly  singular.  This  additional  restriction  on  u(x)  does  not  pose  problems  within 
the  context  of  IgA. 


As  usual,  in  IgA  all  integrals  are  evaluated  in  the  parametric  domain  T.  Thus  numerical  integration  schemes 
have  to  be  developed,  so  that  they  properly  take  into  account  the  geometric  parametrization.  In  particular, 
the  weak  singularity,  established  in  the  physical  domain  T,  should  be  preserved  for  integrals  defined  over 
the  parametric  domain  T.  Further,  since  T  is  partitioned  into  Bezier  elements,  defined  such  that  within 
each  element  the  supported  T-splines  are  C°° ,  locally,  the  map  ip  and  the  basis  functions  are  also  C°° .  This 
is  sufficient  for  developing  numerical  integration  schemes  with  spectral  accuracy  [52].  If  the  integrations 
were  carried  out  in  the  physical  domain,  then  it  would  be  natural  to  use  curvilinear  coordinates  attached 
to  the  surface.  Accordingly,  if  one  uses  parametric  coordinates,  one  needs  to  rely  on  differential  geometry. 
Alternatively,  one  can  avoid  explicit  use  of  curvilinear  coordinates  by  constructing  a  transformation  based 
on  the  singular- value  decomposition  theorem  [19].  This  transformation  is  constructed  as  follows: 

1.  At  a  given  collocation  point,  compute  the  vectors  r1  and  r2  as 

i  _  ^(771,772)  2  _  5^(771,772) 

t  ana  t  ^ 

orji  dr]2 

where  (771,772)  are  the  coordinates  of  the  collocation  point  in  F. 
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2.  Form  a  Jacobian  matrix 


where  the  subscripts  refer  to  the  Cartesian  components  in  a  global  coordinate  system,  in  which  the 
control  points  are  prescribed. 

3.  Compute  the  reduced  singular-value  decomposition 

J  =  UY.V*  , 


so  that  17  is  a  3  x  2  matrix  and  E  is  a  2  x  2  diagonal  matrix  with  no  zeros  on  the  diagonal. 

4.  Define  the  reparametrization  <^(^1,^2)  °f  ^(£17  £2)  via 


^(£i>6)  =  =  <p(V"£  10  ) 


where  £  =  (£i,6)T- 

5.  Compute  the  orthonormal  vectors  f 1  and  f2  as 

=  md  74, 

where  (^1,^2)  are  the  transformed  parametric  coordinates  of  the  collocation  point. 

This  reparametrization  scheme,  leading  to  orthonormal  basis  vectors,  gives  rise  to  optimal  numerical  integra¬ 
tion  schemes  for  weakly-singular  integrals  [52],  [19].  Further,  based  on  numerical  results  presented  in  Section 
6,  it  appears  that  the  reparametrization  scheme  is  critical  for  numerical  integration  for  CAD  parametrizations 
involving  collapsed  edges. 


Numerical  integration  in  the  parametric  domain  is  carried  out  Bezier-element-wise.  If  the  collocation  point 
is  inside  the  element,  then  we  apply  the  reparametrization  scheme  (£1 , £2)  (^17^2)7  which  transforms 

a  rectangular  element  into  a  parallelogram,  on  which  one  can  implement  a  standard  singular  integration 
scheme  using  a  polar  coordinate  transformation  [52].  If  the  collocation  point  is  outside  the  element,  then 
the  integrand  is  regular,  and  one  can  use  Gaussian  quadratures,  provided  that  the  collocation  point  is  not 
too  close  to  the  element.  Typically,  Gaussian  quadratures  work  well  for  points  located  at  a  distance  d  >  3 h, 
where  h  is  an  element  size,  and  both  dimensions  are  defined  in  the  physical  space  T.  If  d  <  3 h  one  can 
carry  out  recursive  element  subdivision,  which  reduces  h  until  d  >  3 h.  Thus,  in  effect,  Gaussian  quadratures 
are  applied  whenever  the  point  is  outside  of  the  element,  but  for  cases  when  d  <  3h  one  needs  to  combine 
Gaussian  quadratures  with  recursive  element  subdivision. 


6  Numerical  Examples 

6.1  Overview 

In  this  section,  we  present  numerical  examples  emphasizing  various  important  mathematical  and  computa¬ 
tional  aspects  of  IgA  of  BIEs.  All  examples  involve  three  shapes:  a  torus,  a  sphere,  and  a  cube.  The  torus 
is  a  C 00 -surface  which  allows  a  C-parametrization;  actually  one  can  show  that  the  torus  parametrization  is 
C°°  but  this  additional  smoothness  is  not  exploited.  The  sphere  is  a  C,oc -surface  with  a  (72-parametrization, 
because  it  involves  collapsed  edges  at  the  poles.  The  cube  is  a  <72-surface  with  a  (72-parametrization.  All 
shapes  were  constructed  using  standard  parametrizations  based  on  16  (torus),  8  (sphere),  and  6  (cube) 
NURBS  patches.  For  each  shape,  the  patches  were  used  to  generate  five  meshes  via  uniform  refinement 
in  the  parametric  domain,  so  that  at  each  level  of  refinement,  each  element  was  divided  into  four.  Figure 
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4  shows  the  two  coarsest  meshes  for  each  shape.  Continuous  basis  functions  of  degree  p  were  constructed 
so  that,  upon  refinement,  they  remained  C locally  and  continuous  over  patch  boundaries.  In  contrast 
to  continuous  basis  functions,  their  discontinuous  counterparts  were  C locally  and  discontinuous  over 
patch  boundaries.  Unless  otherwise  noted,  all  regular  approximations  involved  p  =  2,  and  degree  elevated 
approximations  involved  p  =  3.  The  number  of  integration  points  in  numerical  integration  schemes  involving 
the  polar  coordinate  transformation  was  dependent  on  the  refinement  level.  At  the  coarse  level,  we  used 

5  points  in  the  radial  direction  and  10  points  in  the  angular  direction.  With  each  refinement,  the  number 
of  points  in  the  radial  direction  stayed  the  same,  while  the  number  of  points  in  the  angular  direction  was 
increased  by  three.  For  regular  integrals,  we  used  the  5x5  Gaussian  quadrature  rule  on  each  subelement. 
These  rules  were  established  empirically  and  no  attempts  were  made  to  optimize  them. 


The  majority  of  numerical  examples  involved  manufactured  exact  solutions  in  the  form 


u(x) 


1 

\x  -  x0\  ' 


(27) 


For  each  shape  the  source  point  £o  was  chosen  far  away  from  the  shape  center  (Fig.  3),  so  that  u(x) 
was  an  analytic  function  which  did  not  involve  near-singular  behavior.  These  manufactured  solutions  were 
chosen  in  order  to  demonstrate  the  necessity  of  numerical  schemes  even  for  problems  with  smooth  solutions. 
The  function  u(x)  was  used  to  construct  the  boundary  data  go  and/or  g jy  for  various  BVPs.  After  that, 
appropriate  BIEs  were  solved  numerically  to  reconstruct  the  full  Cauchy  data. 


The  quality  of  numerical  solutions  was  measured  using  the  L2(F)-error  for  the  Cauchy  data.  To  evaluate  the 
order  of  convergence,  we  defined  the  mesh  size  h  as  the  square  root  of  the  area  of  the  largest  Bezier  element 
in  the  physical  space.  The  estimated  order  of  convergence  (eoc)  for  each  refinement  was  computed  as 


eoc  = 


log(g) 


where  the  subscripts  /  and  c  refer  to  the  fine  and  coarse  meshes,  respectively.  For  p  =  2  the  optimal  order 
of  convergence  for  the  L2(r)-error  for  the  Cauchy  data  is  equal  to  p  +  1  =  3. 


Unless  stated  otherwise,  all  arising  algebraic  problems  were  solved  using  a  preconditioned  GMRES  method 
[47]  with  a  tolerance  of  10~12.  Each  preconditioner  was  constructed  as  the  inverse  of  the  interpolation 
matrix  corresponding  to  the  basis  functions.  As  a  result,  we  were  able  to  reveal  the  spectral  properties  of 
the  collocated  operators  and  significantly  reduce  the  iteration  count. 


In  the  remainder  of  this  section,  we  present  six  case  studies.  Each  study  demonstrates  the  importance  of 
a  particular  aspect.  Those  studies  focus  on  (i)  the  recursive  subdivision  scheme  for  near-singular  integra¬ 
tion  (Section  6.2),  (ii)  local  surface  reparametrization  (Section  6.3),  (iii)  exponential  convergence  of  the 
adopted  integration  scheme  (Section  6.4),  (iv)  spectral  properties  of  the  collocated  operators  (Section  6.5), 
(v)  discontinuous  basis  functions  (Section  6.6),  and  (vi)  approximations  for  mixed  BVPs  (Section  6.7). 

6.2  Recursive  subdivision  for  near-singular  integration 

The  objective  of  this  section  is  to  demonstrate  that  the  recursive  subdivision  scheme  (Section  5.4)  for 
evaluating  near-singular  integrals  is  essential.  To  this  end,  we  considered  the  manufactured  pure  Neumann 
BVP  on  the  torus  and  established  that  the  optimal  order  of  convergence  and  accurate  results  can  be  attained 
only  if  the  subdivision  scheme  was  employed.  The  approximate  solutions  to  this  problem  were  obtained  by 
solving  (18)  and  (19),  using  continuous  basis  functions. 
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(a) 


(b) 


T 


xo  =  (0;  0;  60) 


T 


xo  =  (0;  0;  20) 


z 


(c) 


1 


Xo  =  (0.5;  0.5;  10) 


Figure  3:  Three  representative  shapes:  (a)  torus  (inner  radius  r  =  1  and  outer  radius  R  =  3),  (b)  sphere 
(radius  r  =  1),  and  (c)  cube  (edge  length  a  =  1).  The  source  points  for  the  manufactured  solutions: 
xo  =  (0,  0,  60)  for  the  torus,  Xo  =  (0, 0,  20)  for  the  sphere,  and  xq  =  (1/2, 1/2, 10)  for  the  cube. 
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Figure  5  presents  the  A2(T)-error  of  iih  for  two  numerical  integration  schemes,  with  and  without  recursive 
subdivision.  It  is  clear  that  recursive  subdivision  is  necessary  for  attaining  the  optimal  order  of  convergence. 
Furthermore,  recursive  subdivision  significantly  reduced  the  magnitude  of  the  errors.  The  numerical  example 
is  representative  of  the  other  shapes  and  boundary  conditions. 


h 

With 

subdivision 
L2-error  eoc 

Without 
subdivision 
L2-error  eoc 

2.87E-01 
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Figure  5:  T2(T)-errors  for  two  near-singular  integration  schemes:  with  and  without  subdivision. 


6.3  Surface  reparametrization  scheme 

The  objective  of  this  section  is  to  demonstrate  the  importance  of  the  surface  reparametrization  scheme 
(Section  4)  for  evaluating  singular  integrals,  particularly  when  collapsed  edges  are  involved.  To  this  end,  we 
considered  the  manufactured  pure  Dirichlet  BVP  on  the  sphere.  The  approximate  solutions  to  this  problem 
were  first  obtained  by  using  the  SBIE  and  discontinuous  basis  functions.  The  corresponding  linear  algebraic 
problem  is 


Vt  =  ls.  (28) 

Note  that  the  SBIE  requires  one  to  invert  the  matrix  corresponding  to  the  single-layer  operator,  which  is 
not  optimal,  as  far  as  the  spectral  properties  are  concerned.  Nevertheless,  it  allows  us  to  demonstrate  that 
the  surface  reparametrization  scheme  is  necessary  even  for  a  naturally  weakly-singular  operator. 


Figure  6  presents  the  L2(T)-error  for  the  solutions  of  (28),  using  two  numerical  integration  schemes,  with 
and  without  the  reparametrization.  It  is  clear  that  the  reparametrization  is  necessary  for  attaining  the 
optimal  order  of  convergence  and  small  errors.  This  example  is  representative  of  other  CAD  parametrizations 
involving  collapsed  edges. 


Alternatively,  one  can  solve  the  manufactured  problem  using  the  HSBIE  (20)  and  the  discontinuous  basis 
functions.  In  this  case,  similar  to  equation  (28),  the  reparametrization  scheme  is  necessary  for  accurate 
integration  at  collocation  points  near  collapsed  edges.  Further  the  reparametrization  scheme  is  natural  for 
computing  the  tangent  gradient  required  for  regularizing  the  hyper-singular  operator.  Figure  7  presents  the 
L2(r)-error  of  on  the  sphere  using  SBIE  and  HSBIE;  results  for  the  SBIE  are  identical  to  those  presented 
in  Figure  6.  It  is  clear  that  the  two  approaches  yield  similar  results,  and  therefore  both  are  acceptable. 
Thus  we  can  conclude  that  the  surface  reparametrization  scheme  is  necessary,  and  it  is  capable  of  delivering 
optimal  and  accurate  numerical  solutions  even  if  the  hyper-singular  operator  is  involved. 
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Figure  6:  The  manufactured  pure  Dirichlet  BVP  for  the  sphere:  L2(r)-errors  for  two  singular  integration 
schemes,  with  and  without  reparametrization. 
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6.4  Exponential  convergence  of  the  integration  scheme 

The  objective  of  this  section  is  to  demonstrate  that  the  adopted  numerical  integration  scheme  is  exponentially 
convergent  with  respect  to  the  number  of  integration  points.  The  demonstration  involves  all  three  shapes. 
Thus  the  integration  scheme  was  tested  on  problems  involving  non-smooth  surfaces  and  parametrizations. 
For  each  shape,  the  test  problem  was  a  pure  Neumann  BVP  with  the  exact  solution 

u(x)  =  Xi  +  X2  +  X3,  x  £  fl . 

For  this  choice,  one  can  prove  that  the  approximation  error  is  exactly  equal  to  zero  for  every  shape,  and 
therefore  the  chosen  test  problems  are  ideal  for  assessing  numerical  integration  errors. 


Figure  8  presents  the  L2(F)-error  of  the  solution  Uh ■  The  results  were  obtained  by  using  the  coarsest  meshes, 
while  the  number  of  integration  points  in  each  direction  was  increased.  The  results  confirm  an  exponential 
convergence  for  every  shape.  Note  that  we  were  able  to  reach  the  machine  precision  for  the  sphere  and  the 
cube,  while  for  the  torus  the  error  stagnated  near  10“  10  due  to  round-off  errors  in  the  adopted  numerical 
integration  scheme.  If  desired,  this  issue  can  be  resolved  by  using  a  more  sophisticated  singular  integration 
scheme  proposed  in  [19]. 

6.5  Spectral  properties  of  integral  operators 

It  is  well-known  that,  upon  discretization,  the  operators  +  K  and  \X  —  K'  give  rise  to  well-conditioned 
matrices.  The  objective  of  this  section  is  to  demonstrate  that  one  should  choose  the  governing  BIEs  so  that 
one  takes  advantage  of  this  property.  Accordingly,  for  pure  Neumann  BVPs,  the  SBIE  is  a  natural  choice. 
In  contrast,  for  pure  Dirichlet  BVPs,  one  should  choose  the  HSBIE.  This  choice  would  be  problematic 
for  conventional  collocation  BEMs  but  it  is  legitimate  for  IgA.  This  point  will  be  supported  by  numerical 
examples  presented  in  Section  6.5.1.  Further,  we  show  that  for  mixed  boundary-value  problems  one  should 
choose  the  SBIE  on  Tn  and  the  HSBIE  on  r^.  Numerical  examples  presented  in  Section  6.5.2  show  that, 
as  far  as  spectral  properties  and  iteration  counts  are  concerned,  this  choice  is  superior  than  uniform  use  of 
SBIE.  Numerical  examples  involved  the  manufactured  solutions  for  all  three  shapes  and  results  are  presented 
for  both  iteration  counts  and  condition  numbers. 

6.5.1  Dirichlet  problems 

In  principle,  a  pure  Dirichlet  BVP  can  be  solved  using  either  the  SBIE  or  HSBIE.  The  corresponding  linear 
algebraic  problems  are 

Vi-  =  ls 

and 

:')*  =  /„■ 

Thus  the  SBIE-based  approach  requires  one  to  invert  V ,  whereas  the  HSBIE-based  approach  requires  one 
to  invert  \In  —  K' . 

Figure  9  presents  the  iteration  counts  and  spectral  condition  numbers  k  as  functions  of  h  for  the  three 
shapes.  It  is  clear  that  the  HSBIE  is  a  better  choice  than  the  SBIE  both  in  terms  of  the  iteration  counts  and 
spectral  properties.  For  the  torus  and  sphere,  the  iteration  counts  and  spectral  properties  for  the  HSBIE 
are  independent  of  the  mesh  size.  This  is  in  agreement  with  theoretical  results  based  on  compactness  of  K! . 
In  this  regard,  let  us  observe  that  upon  discretization  the  torus  remains  a  C2-surface,  whereas  the  sphere 
becomes  a  <72-surface.  Thus  the  numerical  results  for  the  torus  are  fully  expected,  while  the  results  for  the 
sphere  need  additional  theoretical  considerations.  For  the  cube,  both  the  iteration  count  and  the  spectral 
condition  number  show  a  logarithmic  dependence  on  h.  For  the  cube,  the  mathematical  foundations  are  not 
well-established  because  1C  is  not  compact.  For  the  SBIE,  the  iteration  count  should  grow  as  l/Vh  whereas 
k  should  grow  as  1/h.  Surprisingly,  these  scalings  hold  only  for  the  cube.  It  is  unclear  to  us  why  the  results 
for  the  torus  are  better  than  expected  and  for  the  sphere  worse  than  expected. 
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Numerical  integration  L2(r)-errors  for  the  (a)  torus  (b)  sphere,  and  (c)  cube. 
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Figure  9:  Iteration  counts  of  the  preconditioned  GMRES  method  and  condition  numbers  k  for  the  (a)  torus 
(b)  sphere,  and  (c)  cube  . 
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6.5.2  Mixed  boundary-value  problem 


In  this  section,  we  present  numerical  examples  suggesting  that  one  should  choose  the  SBIE  on  Tjv  and 
the  HSBIE  on  as  opposed  to  using  the  SBIE  on  the  entire  T.  For  our  purposes,  we  choose  boundary 
conditions  as  shown  in  Figure  10.  For  the  torus  and  sphere,  the  Dirichlet  (Neumann)  boundary  conditions 
are  prescribed  on  the  upper  (lower)  halves.  For  the  cube,  the  Neumann  boundary  conditions  are  prescribed 
on  the  top  and  bottom  faces,  and  the  Dirichlet  boundary  conditions  on  the  other  faces.  Similar  to  pure 


I  Dirichlet  BC 
I  Neumann  BC 


Figure  10:  Mixed  boundary  conditions  for  the  torus,  sphere,  and  cube. 

Dirichlet  BVPs,  mixed  BVPs  can  be  solved  with  our  without  the  HSBIE.  We  refer  to  the  former  method  as 
SBIE/HSBIE  and  to  the  latter  one  SBIE/SBIE.  In  the  SBIE/HSBIE,  the  natural  domain  for  the  HSBIE  is 
r d,  as  in  Dirichlet  BVPs,  while  the  SBIE  is  natural  for  Pjv.  In  the  SBIE/SBIE,  the  SBIE  must  be  applied 
in  its  unnatural  domain  Tp.  In  terms  of  linear  algebra,  SBIE/HSBIE  allows  one  to  avoid  inverting  matrices 
with  unfavorable  spectral  properties. 


Figure  11  presents  the  iteration  counts  for  the  mixed  BVPs  for  the  three  shapes.  It  is  clear  that  (i)  the 
SBIE/HSBIE  is  superior  to  SBIE/SBIE,  and  (ii)  the  iteration  counts  grow  with  mesh-refinement  for  all 
three  problems;  for  the  SBIE/HSBIE  scheme,  the  iteration  counts  grow  logarithmically  with  h-1.  We  do  not 
present  comparisons  for  the  spectral  properties  because  such  comparisons  strongly  depend  on  the  definition  of 
spectral  properties.  That  is,  the  spectral  properties  for  the  entire  matrices  are  different  from  those  obtained 
using  Schur  complements;  in  contrast,  the  use  of  Schur  complements  has  minimal  effects  on  iteration  counts. 

6.6  Discontinuous  basis  functions 

The  objective  of  this  section  is  to  demonstrate  that  discontinuous  basis  functions  are  critical  for  approxi¬ 
mating  t  when  it  is  discontinuous.  Those  problems  include  not  only  non-smooth  domains  like  a  cube,  but 
also  mixed  BVPs  defined  on  smooth  domains;  in  the  latter  case,  t  may  be  discontinuous  at  the  interface 
between  T jj  and  Tv-  For  demonstration  purposes,  we  solved  the  manufactured  pure  Dirichlet  BVP  on  the 
cube.  In  this  problem,  t  is  discontinuous  at  the  edges  and  vertices  because  of  the  discontinuous  normal. 
Further,  the  normal  discontinuity  does  not  allow  us  to  collocate  the  operators  —  JC  and  V  at  the  edges 
and  vertices.  This  creates  two  options:  (i)  one  can  use  the  SBIE  with  either  continuous  or  discontinuous 
basis  functions,  or  (ii)  one  can  use  either  the  SBIE  or  HSBIE  with  discontinuous  basis  functions  because  they 
require  collocation  points  off  the  edges  and  vertices.  We  pursue  the  first  option  as  it  allows  us  to  compare 
continuous  versus  discontinuous  basis  functions. 
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Figure  11:  Iteration  counts  of  the  preconditioned  GMRES  method  for  the  (a)  torus,  (b)  sphere,  and  (c) 
cube. 
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Figure  12  presents  the  L2(r)-error  of  th  for  two  approximations,  one  involves  only  continuous  basis  functions 
and  the  other  only  discontinuous  ones.  It  is  clear  that  for  the  continuous  basis  functions  th  converges  very 
slowly.  In  contrast,  for  the  discontinuous  basis  functions,  th  converges  at  the  optimal  rate  and  delivers  very 
accurate  solutions. 
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Figure  12:  L2(r)-errors  for  two  approximations,  with  and  without  the  discontinuous  basis  functions. 


6.7  Approximations  for  mixed  boundary-value  problems 

In  approximation  theory  for  BIEs,  it  is  well  established  that  for  mixed  boundary-value  problems  one  should 
use  different  approximations  for  the  Cauchy  data  u  and  t.  In  particular,  to  attain  the  optimal  convergence 
rate  for  the  Cauchy  data,  approximations  for  u  should  be  one  degree  higher  than  those  for  t  [58].  In  this 
section,  we  present  numerical  results  supporting  this  statement.  Further,  we  present  results  suggesting  that 
approximations  of  the  same  degree  are  capable  of  delivering  the  optimal  convergence  rate  for  the  torus  and 
sphere. 


For  demonstration  purposes,  we  considered  the  same  mixed  BVPs  as  in  Section  6.5.  These  problems  were 
solved  using  the  SBIE/HSBIE  scheme,  discontinuous  basis  functions  for  approximating  t ,  and  continuous 
regular  ( p  =  2)  and  degree  elevated  ( p  =  3)  basis  functions  for  approximating  u.  Figure  13  presents  the 
L2(r)-errors  obtained  using  regular  basis  functions  for  u.  It  is  clear  that  the  L2(r)-error  for  u  converges 
optimally  for  all  three  cases.  In  contrast,  the  rate  of  convergence  for  t  is  optimal  for  the  torus  and  sphere, 
but  not  for  the  cube.  Perhaps  the  results  for  the  torus  and  sphere  are  more  surprising  than  those  for  the 
cube,  as  we  expected  suboptimal  convergence  rates  for  t  for  all  three  cases. 


Figure  14  presents  the  L2(r)-errors  obtained  using  degree  elevated  basis  functions  for  u.  It  is  clear  that  in  all 
three  cases  the  L2(r)-errors  for  both  u  and  t  exhibit  optimal  convergence  rates.  Note  that  for  the  cube  the 
errors  corresponding  to  the  degree  elevated  basis  functions  are  significantly  smaller  than  those  corresponding 
to  the  regular  basis  functions.  For  the  torus  and  sphere,  the  errors  in  u  corresponding  to  the  degree  elevated 
basis  functions  are  also  significantly  smaller  than  those  corresponding  to  the  regular  basis  functions,  which 
is  not  surprising  simply  because  we  used  higher-order  approximations.  In  contrast,  for  the  torus  and  sphere, 
the  degree  elevated  approximation  for  u  had  a  minor  impact  on  the  errors  for  t. 


7  Summary 

In  this  paper,  we  adopted  IgA  as  the  foundation  for  solving  BIEs  corresponding  to  Laplace’s  equation. 
Accordingly,  we  focused  on  problems  defined  on  C2-surfaces,  which  are  common  in  IgA.  Our  theoretical 
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Figure  13:  L2(P)-errors  for  the  flux  t  and  potential  u  on  the  (a)  torus  (b)  sphere,  and  (c)  cube.  The 
approximations  involved  continuous  basis  functions  degree  p  =  2  for  u  and  discontinuous  basis  functions  of 
degree  p  =  2  for  t. 
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Figure  14:  L2(r)-errors  for  the  flux  t  and  potential  u  on  the  (a)  torus  (b)  sphere,  and  (c)  cube.  The 
approximations  involved  continuous  basis  functions  degree  p  =  3  for  u  and  discontinuous  basis  functions  of 
degree  p  =  2  for  t  (solid  lines).  The  dashed  lines  correspond  to  the  results  presented  in  Figure  13. 
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results  and  numerical  schemes  take  full  advantage  of  this  geometric  smoothness.  In  this  regard,  our  work  is 
a  major  departure  from  previous  papers  concerned  with  applying  IgA  to  BIEs,  where  IgA  was  considered 
merely  as  a  BEM  with  different  basis  functions. 


Numerical  schemes  developed  in  this  work  allow  one  to  apply  collocation  schemes  to  both  SBIE  and  HS- 
BIE;  ordinarily  the  HSBIE  can  be  analyzed  within  Galerkin’s  setting  only.  The  access  to  both  SBIE  and 
HSBIE  was  exploited  for  constructing  governing  linear  algebraic  equations  with  optimal  spectral  proper¬ 
ties.  Furthermore,  all  integral  operators  involved  in  the  SBIE  and  HSBIE  were  reduced  to  weakly-singular 
integrals,  for  which  we  adopted  effective  numerical  integration  schemes  available  in  the  literature,  includ¬ 
ing  the  polar  coordinates  transformation,  local  reparametrization  of  the  surface,  and  recursive  subdivision. 
This  combination  resulted  in  an  exponentially  convergent  integration  scheme,  essential  for  attaining  optimal 
approximation  properties. 


Adopting  IgA  for  numerical  treatment  of  BIEs  necessitated  the  introduction  of  the  following  numerical 
schemes: 

•  Discontinuous  basis  functions  for  treating  edges  and  vertices  of  non-smooth  domains,  and  irregularities 
associated  with  collapsed  edges  and  extraordinary  points,  both  common  in  IgA. 

•  Local  reparametrization  of  surfaces  necessary  for  effective  and  accurate  integration  schemes,  especially 
in  the  vicinity  of  collapsed  edges. 

•  Degree  elevated  basis  functions  advantageous  for  solving  mixed  boundary-value  problems. 


There  are  several  open  issues,  whose  resolution  may  significantly  advance  IgA  of  BIEs.  Among  them  are: 

•  Mathematical  foundations  for  collocation  schemes  for  BIEs  corresponding  to  pure  BVPs  defined  on 
non-smooth  domains. 

•  Mathematical  foundations  for  collocation  schemes  for  BIEs  corresponding  to  mixed  BVPs  defined  on 
smooth  and  ultimately  non-smooth  domains. 

•  Efficient  integration  schemes  for  Galerkin  discretizations. 

•  Better  numerical  integration  schemes  that  exploit  smoothness  of  the  basis  functions. 

•  Adaptive  numerical  integration  and  approximation  schemes,  in  the  spirit  of  h  —  p  —  k  methods,  intro¬ 
duced  and  described  herein  [24,  25]. 

In  conclusion,  let  us  mention  that  most  of  the  theoretical  results  and  computational  numerical  schemes  can 
be  extended  to  equations  of  linear  elasticity  and  Stokes  equations  of  fluid  mechanics.  Those  extensions  may 
have  a  significant  impact  on  the  development  of  tools  for  engineering  design  and  analysis. 
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A  Compactness  of  double-layer  operators  on  C2-surfaces 

In  this  section  we  prove  that  the  double-layer  operators  K,  and  1C  defined  on  T  G  C2  are  compact  operators 
on  the  space  of  continuous  functions. 


Consider  a  bounded  domain  U  C  R3 4  with  P  :=  dO.  By  definition  r  G  C2  if  it  satisfies  the  following 
conditions: 

1.  For  any  x  GT  there  exists  a  constant  R  >  0,  independent  of  x,  and  a  neighborhood  Nx  C  R3  such  that 
dist(a;,  T\NX)  >  R- 

2.  The  surface  :=  P  (~l  Nx  is  an  image  of  a  domain  Nx  C  R2  under  the  map  f)x, 

3.  The  map  ifx  is  bijective  and  twice  continuously  differentiable; 

4.  The  maps  ipx  and  fjf1  are  Lipschitz  continuous. 

Theorem  A.l.  IfTG  C2  then  there  exists  a  constant  C  >  0  such  that 

(x-y)-n(y)  <  C 
4tt\x  —  y  |3  —  \x  —  y\ 

for  all  x,y  GT. 
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Proof.  Let  x  &  T.  If  'y  f  Tx,  then  there  exists  a  constant  R  >  0  such  that  dist(x,  r\rx)  >  R  and  consequently 
| y  —  x\  >  R.  As  a  result  we  obtain  the  estimate 

(x  —  y)  ■  n(y)  diam(fl)  diam(fl)2  1 
47r|a;  —  y  |3  —  47t_R3  —  47tI?3  \x  —  y\ 


If  y  €TX  then  there  exist  x  €  Nx  and  y  £  Nx  such  that  ipx(x)  =  x  and  f’xiv)  =  V  because  ipx  is  bijective. 
Since  tpx  is  twice  continuously  differentiable,  we  can  use  the  following  expressions  for  the  normal 


n(y) 


&MS)  v  di/ixjy) 
dyi _ dy2 

d^x(y)  v  ^bx(y) 
dy  1  dy2 


and  Taylor’s  expansion 

x-y  =  tpx{x)  -  ipx(y)  =  ^w^-ix  i  -  yi)  +  (x2  -  y2)  +  0( \x  -  y\2) 

oyi  dy2 

to  obtain  the  estimate 


(x-y)-n{y)  =  OQx  -  y\2). 

Since  iff1  is  Lipschitz  continuous,  there  exists  a  constant  l  >  0  such  that  l\x  —  y\  <  \ipx(x)  —  il>x(y)\  =  \x  —  y\, 
which  implies 


{x  -  y)  ■  n{y)  =  Q  f  1  \ 

47r|a;-j/|3  \\x-y\J 


□ 


This  theorem  implies  that  the  operator  K  is  weakly  singular.  Similarly,  we  can  establish  the  bound 


(x-y)-  n(x) 
47r|a:  —  y  |3 


<  |  °  |  Vx,y  G  T 

\x  -  y\ 


and  thus  establish  that  the  operator  fC  is  also  weakly  singular. 


Using  well-established  techniques,  it  can  be  proved  that  the  operators  /C,  /C'  :  C(T)  ->  C'(T)  are  compact  and 
the  operators  \l  +  K,  and  \l  —  1C  are  invertible  [30,  Chapter  3].  Furthermore,  if  we  assume  that  there  exists 
a  bounded  interpolation  operator  that  maps  continuous  functions  onto  continuous  basis  functions  (24),  one 
can  invoke  compact  perturbations  theory  to  prove  that  the  discrete  problems  (18)  and  (20)  are  uniquely 
solvable  and  their  solutions  iq,  and  th  satisfy  quasi-optimal  error  estimates  in  the  L°°(r)-norm  [18,  Chapter 
XII.  1], 
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